% @author:夜剑听雨
% E-mail:2530595378@qq.com
% 2025.08.29
% 绘制反演台阶速度图
clc; clear; close all;
load('./save_dir/Inv/PSO_Inv_model3.mat')
PSO_Inv_model = Inv_model;

%% 台阶速度图
% 模型1：速度递增，模型2：低速夹层，模型3：高速夹层
% Vs_true = [150 300 400 600]; h_true = [2 4 6];
% Vs_true = [250 120 350 500]; h_true = [4 2 4];
Vs_true = [100 200 500 400]; h_true = [4 2 4];
True_mod = [Vs_true h_true];

PSO_mean = mean(PSO_Inv_model);
mean_err = mean((abs(PSO_mean-True_mod)./True_mod)) * 100;

%% 画图
maxDepth = 20;
layNum = 4;
% 反演初始模型台阶速度图
figure;

% 1.绘制上下边界（不添加图例）
% 模型1上下边界
% Vs_ub = [300 400 500 800];h_ub = [3 5 8];
% Vs_lb = [100 200 300 500];h_lb = [1 2 4];

% 模型2上下边界
% Vs_ub = [300 300 500 600];h_ub = [6 3 6];
% Vs_lb = [100 100 300 400];h_lb = [2 1 2];
% 模型3上下边界
Vs_ub = [150 300 600 500];h_ub = [6 3 6];
Vs_lb = [50 100 400 300];h_lb = [2 1 2];

[vsLb, depthLb] = calcStairsData(Vs_lb, h_ub, maxDepth);
[vsUb, depthUb] = calcStairsData(Vs_ub, h_lb, maxDepth);

% 绘制边界线（灰色虚线）
stairs(vsLb, depthLb, 'Color', [0.5 0.5 0.5], 'LineStyle', '--', 'LineWidth', 1); hold on
stairs(vsUb, depthUb, 'Color', [0.5 0.5 0.5], 'LineStyle', '--', 'LineWidth', 1); hold on

% 2.绘制真实值（黑色实线，不添加图例）
[vsTrue, depthTrue] = calcStairsData(Vs_true, h_true, maxDepth);
stairs(vsTrue, depthTrue, 'k-', 'LineWidth', 1.5); hold on

%% ====== 为反演结果分配不同颜色 ======
% 创建5种不同颜色（使用lines颜色图）
colors = lines(5);  % 生成5种鲜明对比的颜色

% 初始化图例句柄和文本
legendHandles = [];  % 句柄向量
legendTexts = {};    % 标签元胞数组

% 3. 绘制PSO反演结果并添加到图例
for i = 1:5
    PSO = PSO_Inv_model(i,:);
    err = mean(abs(PSO - True_mod) ./ True_mod) * 100;
    [vsInv, depthInv] = calcStairsData(PSO(1:layNum), PSO(layNum+1:end), maxDepth);

    % 使用不同的颜色绘制每个反演结果
    p_inv = stairs(vsInv, depthInv, '--', 'Color', colors(i,:), 'LineWidth', 1);
    hold on

    % 将反演结果添加到图例
    legendHandles = [legendHandles, p_inv];
    legendTexts{end+1} = sprintf('Inv%d: %.2f%%', i, err);
end

% 创建图例（只包含5个反演结果）
lgd = legend(legendHandles, legendTexts, ...
    'Location', 'southwest', 'FontName', 'Microsoft YaHei', 'FontSize', 8);
lgd.Color = 'none';

% 设置y轴的范围和刻度
ylim([0 maxDepth]); yticks(0:5:maxDepth);
% 设置x轴的范围和刻度
xlim([0 800]); xticks(0:100:800);
% 设置坐标轴线宽
set(gca, 'LineWidth', 1);
% 反转y轴
set(gca, 'XAxisLocation', 'bottom');
set(gca, 'YDir', 'reverse');

title('Inversion model', 'FontSize', 12, 'FontName', 'Microsoft YaHei');
xlabel('S-wave velocity (m/s)', 'FontSize', 12);
ylabel('Depth (m)', 'FontSize', 12);

% 优化布局
box on;
set(gca, 'LooseInset', get(gca, 'TightInset') + 0.01);

%% 画平均值图

figure;

[vsLb, depthLb] = calcStairsData(Vs_lb, h_ub, maxDepth);
[vsUb, depthUb] = calcStairsData(Vs_ub, h_lb, maxDepth);

% 绘制边界线（灰色虚线）
p_lb = stairs(vsLb, depthLb, 'Color', [0.5 0.5 0.5], 'LineStyle', '--', 'LineWidth', 1); hold on
stairs(vsUb, depthUb, 'Color', [0.5 0.5 0.5], 'LineStyle', '--', 'LineWidth', 1); hold on

% 2.绘制真实值（黑色实线，不添加图例）
[vsTrue, depthTrue] = calcStairsData(Vs_true, h_true, maxDepth);
p_true = stairs(vsTrue, depthTrue, 'k-', 'LineWidth', 1.5); hold on

% 3. 绘制PSO反演结果并添加到图例
[vsInv, depthInv] = calcStairsData(PSO_mean(1:layNum), PSO_mean(layNum+1:end), maxDepth);
p_inv1 = stairs(vsInv, depthInv, 'r-', 'LineWidth', 1.5); hold on

% 预分配图例句柄和文本
legendHandles = cell(3, 1); % +1 给真实值
legendTexts = cell(3, 1);
legendHandles{1} = p_lb;
legendTexts{1} = 'Boundary';
legendHandles{2} = p_true;
legendTexts{2} = 'True';
legendHandles{3} = p_inv1;
legendTexts{3} = sprintf('PSO: %.1f%%', mean_err);

% 
lgd = legend([legendHandles{1},legendHandles{2}, legendHandles{3}], legendTexts,...
    'Location', 'southwest', 'FontName', 'Microsoft YaHei', 'FontSize', 8);
% 设置图例为半透明
lgd.Color = 'none';          % 完全透明背景

% 设置y轴的范围和刻度
ylim([0 maxDepth]); yticks(0:5:maxDepth);
% 设置x轴的范围和刻度
xlim([0 800]); xticks(0:100:800);
% 设置坐标轴线宽
set(gca, 'LineWidth', 1);
% 反转y轴
set(gca, 'XAxisLocation', 'bottom');
set(gca, 'YDir', 'reverse');

title('Inversion model', 'FontSize', 12, 'FontName', 'Microsoft YaHei');
xlabel('S-wave velocity (m/s)', 'FontSize', 12);
ylabel('Depth (m)', 'FontSize', 12);

% 优化布局
box on;
set(gca, 'LooseInset', get(gca, 'TightInset') + 0.01);